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In this work, we consider excited many-body mean-field states of bosons in a double-well optical 
lattice by investigating stationary Bloch solutions to the non-linear equations of motion. We show 
that, for any positive interaction strength, a loop structure emerges at the edge of the band structure 
whose existence is entirely due to interactions. This can be contrasted to the case of a conventional 
optical (Bravais) lattice where a loop appears only above a critical repulsive interaction strength. 
Motivated by the possibility of realizing such non-linear Bloch states experimentally, we analyze the 
collective excitations about these non-linear stationary states and thereby establish conditions for 
the system's energetic and dynamical stability. We find that there are regimes that are dynamically 
stable and thus apt to be realized experimentally. 



I. INTRODUCTION 

A wealth of non-trivial quantum states have been re- 
alized with Bose-Einstein condensates (BECs) in optical 
lattices [1, 2]. Conventionally, these states occupy the 
lowest energy configuration of the parent Hamiltonian of 
the system. On the other hand, motivated by the low 
levels of dissipation in ultracold atomic gases, consider- 
able recent effort has shifted to the possibility of realizing 
interesting quantum states not as the ground state, but 
as long-lived excited states [3]. In particular, recent ex- 
periments have realized condensates in excited bands [4- 
8]. Such efforts open the door to achieving, for instance, 
states having Neel order, which is notoriously difficult to 
realize as the ground state of a system of bosons [9, 10], 
and novel fermionic states [11, 12] that have no paral- 
lels in solid state physics. Unconventional excited states 
realized with ultracold atoms have been considered the- 
oretically in Refs. [13-18]. 

In the presence of interactions, the nonlinear Bloch 
band of a BEC can develop interesting features. Al- 
though its ground state is fairly conventional, at the Bril- 
louin zone boundary, the band can develop a cusp and 
subsequently form a loop as interaction is increased [19- 
22]. However, for Bravais lattices (with one lattice site 
per unit cell) such a loop appears only above a critical 
interaction strength, which can be large. For all realistic 
situations, the loop is also small and hence difficult to 
detect. This limits the proposed experimental detection 
of the loops to only indirect signals such as the nonlinear 
Landau-Zener effect [19, 23, 24]. 

Similar looped band structures also appear for BEC on 
a double-well optical lattice [22, 25]. Unlike the conven- 
tional looped structure described above, however, we find 
that a significantly large loop is induced for any interac- 
tion strength, and a large energy separation from the ex- 
cited band is possible by suitably tuning the lattice depth 
and the lattice staggering. With the recent experimental 
realization of double- well optical lattices [5, 26, 27] and 
the concomitant theoretical investigations of the many- 
body phenomena in them [28-35] , it is now an appropri- 
ate moment to consider the possibility of detecting the 



looped band on a double- well lattice. 

Previously, the interaction-induced loop structure 
in double-well optical lattice was considered theoret- 
ically for some specific cases. In Ref. [22], a one- 
dimensional Kronig-Penney potential was used to demon- 
strate period-doubling in a double-well lattice in one di- 
mension. The special form of the potential allowed an- 
alytic solutions to be obtained. A tight-binding model 
of a one-dimensional double-well system was analyzed in 
Ref. [25]. Ref. [36] found an analytical solution at the 
band edge for a specific interaction energy by employing 
the Thomas-Fermi approximation. 

In this work, we compute the interacting Bloch 
solutions of the Gross-Pitaevskii equation for two- 
dimensional double-well lattices using realistic lattice po- 
tentials. This allows us to elucidate the behavior of the 
loop structure at the band edge as a function of the lattice 
potential. Essential to the experimental realization of the 
loop states is their stability. By analyzing the behavior 
of the collective modes about the mean-field solutions, 
we find that a range of states in the excited band are 
dynamical stable, which are in experimentally accessible 
parameter regimes. 

The paper is organized as follows. In Sec. II, we set 
up the problem and indicate the specific lattice used in 
our analysis. We then describe the method of obtaining 
non-linear Bloch solutions to the mean field equations of 
motion, and present the loop band structure in Sec. III. 
In Sec. IV, we compute the collective excitations about 
the non- linear Bloch solutions. These are used to estab- 
lish the energetic and dynamical stability of the system 
prepared in these states. In Sec. V, we present a tight- 
binding model that captures some, but not all, of the 
salient features of the more-accurate continuum results. 
We conclude in Sec. VI. 



II. THE DOUBLE- WELL OPTICAL LATTICE 

In a recent experiment [26], a double- well lattice was 
realized by superimposing standing waves in the x and 
y directions. To ensure phase stability, a single light 
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Figure 1: (Left) The double-well lattice potential used in this 
work. The unit cell (dotted red square) and lattice constant 
a are shown. Each black cross corresponds to a site of a 
bosonic annihilation operator in Eq. (8). The minima are 
colored deep-blue (gray) , in accordance to the scale on lower- 
right. (Upper-right) A slice of the lattice potential along the 
diagonal line on the left panel. 



source was used and the lattice is formed in a folded- 
retroreflected geometry. We consider a lattice potential 
that is realizable with this geometry, given by 



III. MEAN-FIELD ANALYSIS OF THE 
INTERACTING BAND STRUCTURE 

In this Section we describe the numerical method used 
to obtain periodic mean-field stationary states of Eq. (2). 
These solutions are shown to exhibit extra "looped" states 
which are then further analyzed. We concentrate on 
the semi-classical regime of H. and approximate tp with 
if> = (ijj). Thus, we seek solutions of the Gross-Pitaevskii 
equation (GPE) 

Mr) = Jv^(r) + Mat(r)^(r) + g \^(r)\ 2 V(r) (3) 

Among the solutions of this nonlinear differntial equa- 
tion, we are interested in the solutions of the Bloch form: 

V>nk(r) = e ikr u nk (r) (4) 

where n is the band index, k is the crystal momentum, 
and w„k( r ) nas the peridodicity of the underlying lattice. 
The corresponding mean-field energy per unit cell is then 



V\Ax,y) = Vi [cos(2k L x) - cos{2k L y)) 
+2V 2 cos(2k L x) cos(2k L y) 



(1) 



where k^ — ir/a and a — A/%/2 is the lattice constant 
(shown in Fig.l). By adjusting the path- length differ- 
ences of the lattice beams, a more general lattice can be 
generated. We further restrict ourselves to the case where 
2V2 > V\ > which gives double- well lattices with de- 
generate maxima and staggered local minima (see Fig. 1) 

With this lattice potential, the full Hamiltonian of the 
system is given by 
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where g 



ip(r) describes the destruction of a 



boson at position r, m is the mass of the constituent 
bosons, and a s > is the s-wave scattering length. 

For the mass and scattering length, we use parame- 
ters for 87 Rb. We consider spatially averaged densities 
p below 1 x 10 14 cm~ 3 , since for larger densities three- 
body losses become important. In the following analy- 
sis we restrict our attention to the case where V\ and 
V2 are less than IOEr, and set the recoil energy to be 

h 2 k 2 

Er = g-r* = h X 1.75 kHz, given by the experimental 
parameters. For such parameters the tight-binding ap- 
proximation is not necessarily valid. For this range of V\ 
and V2, the ground state is in the superfluid phase, well 
away from the Mott insulating transition. 
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(5) 



Throughout the paper, we choose the zero of energy such 
that the ground state energy is zero. 

We numerically compute the nonlinear Bloch band 
by expanding u„k(x) as u„k(x) = ^K CnKeiK r wnere 
K are reciprocal lattice vectors and the summation in- 
cludes a sufficient number of harmonics to ensure accu- 
racy. This is then substituted to Eq. (3) and a set of 
coupled equations is obtained by equating coefficients of 
equal harmonics. Together with the normalization con- 
dition, {c„k} and /i are then solvcd-for numerically. We 
start with an initial solution {c„k} at k = 0, found with 
imaginary-time propagation. Then k is changed stepwise 
and {c„k} is computed through numerical root-finding at 
each step. The non-linear band structures for several sets 
of parameters are shown in Fig. 2. 

The most prominent feature is the emergence of a 
looped band structure when a lattice staggering V\ ^ 
is introduced. Similar phenomena were also discovered 
in interacting BECs on regular ID [19-21, 37, 38] and 
2D [39, 40] lattices. However, a substantial looped band 
on a regular (unstaggered) lattice requires the interaction 
energy to be much larger than the lattice depth, which is 
experimentally a stringent condition. In contrast, a sub- 
stantial looped band on a staggered lattice only requires 
the interaction energy to be greater than the staggering, 
which is achievable as a staggering of < 1% can be real- 
ized. 

For a qualitative understanding of the emergence of 
the looped band on a double-well lattice, we consider 
the band structures for different staggering in Fig. (2). 
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(a) T X M r 




(b) r x m r 




Figure 2: The lowest two bands computed with: gp = OAEr 
and V% = 4.8-Er, while Vi are: (a) (i.e. no staggering), (b) 
0.04_E_r, and (c) 0.2Er. The labelling of the crystal momenta 
k follows the convention where F = (0, 0), X — (it /a, 0) and 
M = (it /a, it /a). 




Figure 3: The wavefunction and current flow in a unit cell 
(shown in Fig. f) for three eigenstates where (a) ko = 
(0.187T, 0), (b,c) ka = (tt,0). The corresponding points on 
the energy band are shown on the upper-left panel. [The 
other degenerate solution at point (b) of the band is given 
by the horizontal reflection of Fig. b] The other panels show 
the contours of [V'WI an d the field of j. The parameters are: 
gp = OAEr, Vi = 0AE R and V 2 = SE R . 



First consider the case of V\ — (Fig. 2a), in which case 
our unit cell is twice the lattice's natural unit cell. It 
is known that an extra solution exists at the band edge 
[22] and a staggering splits the upper and lower band and 
hence a loop is induced (Fig. 2b). Further increasing the 
staggering enlarges the energy gap between the excited 
band and the loop (Fig. 2c). The importance of using 
a double-well lattice is apparent in Fig. 2c, where we 
see a substantial loop well separated from the excited 
band. Another noticeable difference from a conventional 
loop is that it exhibits quadratic band crossing at ka = 
(tt, 71"). This feature is relevant for the nonlinear Landau- 
Zener effect, which was proposed to detect the loop band 
structure [19, 23, 24] experimentally. 

The superfluid current density is obtained from the 
Bloch solution i/)„k(x) as 

jnk = — Im (V4k v- 0nk) 
m 

= - V (K + k) <K'CnK cos (K' • r - K • r) (6) 

KK' 

In Fig. 3 we plot the wavefunction and current den- 
sity of several indicated states. A currentless state is 
found at the band edge of the interaction-induced state 
(Fig. 3c). For a non-interacting system, it is well known 
that the cell-averaged current satisfies j n k = |Vk-E n k- 
This relation, in fact, can also be shown to hold for the 
nonlinear system in Eq. (3). Therefore currentless so- 



lutions can only appear at the T point or the (c) point 
in Fig. (3), where the energy band attains its local ex- 
trema. The state at (c), however, has identically zero 
current everywhere instead of just an average zero cur- 
rent. This is because the period doubled state can be 
taken to be real everywhere, and hence currentless by 
virtue of Eq. (6). It is easy to see that a real solution is 
possible only at the T, X or M points of the Brillouin 
zone. Because ip(r) — ^ K CKe*' K+k ' r , the reality of ip 
requires 2k be equal to a reciprocal lattice vector and 
that ck = c_K-2k- This is satisfied by the states of the 
loop at the X and M points. 



IV. STABILITY 

A crucial factor determining the realizability of the 
states on the loop is their stability. In general, there are 
two qualitatively different types of instabilities that could 
potentially occur for interacting BEC: energetic instabil- 
ity and dynamical instability. Energetic instability oc- 
curs if the system is not at a local minimum of the mean- 
field energy. This, however, is often unimportant in ex- 
periments since the timescales for energy dissipation are 
long. In contrast , when the system has collective fluctua- 
tions with complex frequencies, small perturbations will 
grow exponentially fast: a dynamical instability. This 
type of instability does not require energy dissipation, 
and it will cause rapid depletion and fragmentation of 
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Figure 4: (Color online) Nonlinear Bloch bands with the 
same paramters of Fig. (3), in the first quadrant of the Bril- 
louin zone. States in blue region are stable while states in 
gray regions are unstable, both dynamically and energetically. 
Red regions are energetically unstable but dynamically stable. 
The excited band is separated from the ground band by about 
IEr and is not shown here. 



the condensate [41, 42]. Thus realizing a dynamically 
unstable state is difficult if the instability timcscale is 
much shorter than the experimental timescale. 

To analyze the stability of states, we expand the Hamil- 
tonian to quadratic order in the field: ip — > rp n ^ + 
e lk r 0(x) where ipnk is the stationary non-linear Bloch 
solution. Hence the term linear in <p vanishes and the 
term of second order in d> is: 



Figure 5: Stability of the eigenstate at ka = (tv, 0) on the 
loop, as a function of V\ (the staggering) and Vi (the lattice 
depth) in Eq. (1), with gp = OAEr (corresponding to a den- 
sity p — 9 x 10 13 cm~ 3 ). The hashed region has no loop. The 
shaded region is dynamically unstable and the white region 
is energetically unstable but dynamically stable. Contours 
show the loop size (i.e. vertical distance between (b) and (c) 
in Fig. 3) in units of Er. 



suppressed; (2) as the staggering becomes smaller, a loop 
is formed but with all states on it dynamically unstable; 
(3) further reducing staggering enlarges the loop and a 
band of stable states is formed on it, of which Fig. 4 is 
an example. 



TIGHT-BINDING MODEL 
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In this Section, we use a tight-binding model to un- 
derstand some of our previous results, in the limit of 
K^q+K' large lattice potentials. The tight-binding Hamiltonian 
we consider is 

(7) 

blb r > + h.c. 



K ' 



-q-K' 



H 



TB 



where q is the Fourier transform of 0(x). p = |V>nk| is 
an effective potential created by the BEC, and p = tp^w 
The collective modes are then numerically found by 
canonical transformation of Eq. (7) . For each crystal mo- 
mentum k, we compute the spectrum w q . The presence 
of any complex u; q implies both dynamical and energetic 
instabilities, while negative w q indicates an energetic in- 
stability. The stabilities of the Bloch states are shown 
in Fig. 4. We find the top of the loop has a region of 
dynamically stable states, which trace out a roughly cir- 
cular path in the Brillouin zone. 

The stability phase diagram of the state at the tail 
at ka = (tt, 0) as a function of V\ and V2 is shown in 
Fig. 5. It is divided into three regions: (1) when the 
lattice staggering is sufficiently large, the loop is totally 



(rr'> 
U 
"2 



■ iQ - r blL 



(8) 



where b r annihilates a boson at position r on a square 
lattice with unit lattice constant, X)( r r') indicates sum- 
mation over nearest neighbors, and Q = (tt, 7r). The 
hopping amplitude is given by t, A > is the on-site en- 
ergy staggering between sites, while U denotes the on-site 
interaction energy between bosons. This model describes 
Eq. (2) in the limit of a strong lattice potential, where 
we associate a bosonic annihilation operator to each local 
minimum of the lattice (see Fig. 1). Note that a term of 
staggered hopping is unnecessary because the lattice has 
4-fold rotational symmetry and all nearest-neighbor links 
are equivalent. 
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In the superfluid regime, the system is described by 
a mean-field equation pb T = —tY^gbr+s + Ae lCi ' r b r + 
U \ b r \ b r , where b r = (b r ) and S denote the positions of 
the nearest-neighbors. This is solved by 

& r = ^(co S ^W Qr sm^)e* r (9) 

where ak satisfies 

2t (cos k x + cos k y ) sinak+A cos a+Upsm «k cos ak = 0. 

(10) 

This yields the nonlinear Bloch band E^/p — Asinak — 
It (cos k x + cos k y ) cos ak — Up (g cos 2 ak — l) • 

In vanishing interaction (U — 0), this tight-binding 
model generates two bands. When U is raised above a 
critical threshold, however, four solutions appear in some 
regions of k near the band edge. These are the extra 
states of the loop. Since Eq. (10) has four solutions if 

(Up) 2/3 > A 2 / 3 + (2t(cosfe x +cosfc y )) 2/3 , the condition 
to have a loop without it filling up the entire Brillouin 
zone (hence forming a separate band) is therefore A < 

Up < ^A 2 /' 3 + (4i) 2/3 ) 3 2 . Further, by demanding all 

Bogoliubov modes to be real, we find that the condition 
of having the state at k = (tz, 0) on the loop dynamically 

stable is t < Since in our lattice [Eq. (1)] V\ 

controls the amount of staggering and V<z controls lattice 
depth, the phase diagram in the continuum (Fig. 5) is in 
qualitative agreement with the tight-binding model. The 
wavefunction at the band edge (cos k x + cos k y = 0) can 
be analytically solved. When Up > A, the state on the 
loop has a — —\. Thus the state at the tail is b r = y/2p 
if e zC *' r = —1 and zero otherwise, i.e. has all particles 
on the lower wells. This too agrees with the continuum 
calculations in Fig. 3. 

This loop structure could also be understood from the 
period-doubled solution [22] . It is known that this tight- 
binding model without staggering (A = 0) admits period- 
p solutions, where p is any positive integer. In particular, 

a period-doubled solution ^ = 4t ^ cos k 2 ^p° s ky ^ +Upal- 
ways exists near the band edge [21]. If we introduce small 
staggering A, the extra band splits in a manner that gen- 
erates a looped band structure (Fig. 2). As the staggering 
gets larger, the splitting between the two bands also in- 
creases but simultaneously suppressing the loop size. As 
splitting becomes large A > Up, the whole loop structure 
is destroyed. 

The adequacy of tight-binding models in describing 
our loop structure could be compared with the situations 
in purely interaction-induced loop structure [20, 40]. In 
those systems, a tight-binding model can never produce 
the loop structure in continuum regardless of the lattice 
depth, and this was attributed to an inappropriate choice 
of Wannier functions [40] . In our system, although a loop 



structure can be captured with the tight-binding model, 
some qualitative features of the continuum calculations 
are missing in the tight-binding model. For instance, the 
tight-binding model predicts that a dynamically stable 
region on the loop, if present, must include the band 
edge. This is not consistent our continuum calculations 
in Sec. Ill, as shown in Fig. 4. 

Although the tight-binding model is simple and gives 
a physical picture, the continuum calculations are more 
relevant from an experimental viewpoint. This is be- 
cause the experimentally tunable parameters are V±, V% 
and p, and it is difficult to compute from these the suit- 
able tight-binding parameters t, A and U of Eq. (8) as 
shown in [43]. In contrast, the results we obtained for 
continuum Hamiltonian Eq. (2) can be directly compared 
with experimental results. There are also experimentally- 
relevant regimes considered above, where a single-orbital 
tight-binding approximation is not correct. 

VI. CONCLUSIONS 

In this work, we have numerically computed the non- 
linear Bloch band structure for an interacting BEC on 
a double- we 11 optical lattice. We also computed the Bo- 
goliubov modes about all the states, and thereby mapped 
out the stability phase diagram as a function of lattice 
depth and staggering for experimentally realistic param- 
eters. A tight-binding model was also considered and 
it is found to reproduce some qualitative features of the 
system. 

We find that the interaction energy required to create a 
looped band is smaller on a double- we 11 lattice, compared 
with a Bravais lattice. Further, we find a stable state on 
the loop can be realized within an experimental accessible 
parameter range. This raises the possiblity of directly 
exciting a BEC onto the loop, or detection of the exotic 
band structure by Raman spectoscopy. 

In future work, a more careful treatment of the im- 
portant possibility of Mott correlations about the states 
on the loop will be interesting to study. When the sys- 
tem is in the tight-binding regime, the Gutziller method 
(which is more general than the treatment provided in 
Sec. V) will be suitable for this purpose. The effects of 
the non-linear Bloch states on the dynamics of BECs in 
double- well lattices will also be of interest. 
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